## Prioritization preferences for COVID-19  -------------------------
## vaccination are consistent across five countries -----------------
## Reproduction materials (Supplementary Appendix) ------------------

set.seed(1234)

# load packages
source("code/helpers.R")

# load data
ranking_df <- readRDS("data/ranking_df.rds")
conjoint_df <- readRDS("data/conjoint_df.rds")
plot_helper_df <- readRDS("data/plot_helper_df.rds")


## -----------------------------------------------------------------
## SF7. Attribute importance (Random forests) ----------------------
## -----------------------------------------------------------------

# random forest estimation function
est_randomforest <- function(df){
  df_model <- df %>%
    select(pat_pref, Age, EarlyRegistration, Gender, Children, Job, Precondition, Citizenship) %>%
    na.omit()
  
  
  X <- model.matrix(pat_pref ~ Age + EarlyRegistration + Gender + Children + Job + Precondition + Citizenship,df_model)
  y <- factor(df_model$pat_pref)
  fit <- randomForest(X[,-1],y)
  var_import <- importance(fit)
  
  df_var_import <- data.frame(var_import)
  df_var_import$var <- rownames(df_var_import)
  rownames(df_var_import) <- NULL
  
  return(df_var_import)
} 

# countrywise estimation
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_randomforest(.))

# dataframe 
figsf7_df <- dplyr::left_join(df_res, plot_helper_df %>% dplyr::distinct(names, attr, levl), by = c("var" = "names")) %>% 
  dplyr::mutate(attr = stringr::str_replace(attr, "Children", "Has children"),
                var_lab = paste(attr, levl, sep=" "),
                var_lab = reorder_within(var_lab,MeanDecreaseGini,country, sep= " "),
                labels_var = paste(attr, levl, sep=" "))

# plot
supp_fig_7 <- ggplot(figsf7_df) + 
  geom_point(aes(x=var_lab,y=MeanDecreaseGini), size=.8) +
  theme_bw() +
  coord_flip() + facet_wrap(~ country, scales = "free") +
  scale_color_grey(name="Respondent PID") +  theme(text = element_text(size=18)) +
  scale_x_discrete(breaks = figsf7_df$var_lab, labels = figsf7_df$labels_var) +
  xlab("") + ylab("Mean Decrease Gini") +
  theme(legend.position = "none", text = element_text(size=14),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        panel.spacing = unit(1, "lines"),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_7

ggsave("figures/supp-fig-7-var-importance.png", plot=supp_fig_7, height=9, width=12)


## -----------------------------------------------------------------
## ST2. Open question analysis -------------------------------------
## -----------------------------------------------------------------

# language specific dfs --------
en <- readLines("data/open-question/EN-orig.txt")
de <- readLines("data/open-question/DE-translated.txt")
br <- readLines("data/open-question/BR-translated.txt")
it <- readLines("data/open-question/IT-translated.txt")
pl <- readLines("data/open-question/PL-translated.txt")

# helper lists ---------
# phrases to tokenize
phrases <- c("has children", "have children", "with children", "children to raise", "had kids" ,"with kids", "has small children",
             "having children", "has younger children", "over 70", "bill gates", "family man", "who contribute", "purchasing power",
             "us born", "over 65", "senior citizen", "senior citizens", "nursing home", "want it", "sign up", "want the vaccine", "not taken", 
             "age of children", "expected gains", "I will never do it", "without vaccines",
             "world leaders", "value in society", "system relevant", "vaccinations are way too premature", "everyone who wants",
             "benefit to the general public", "first come", "function for society", "guinea pig", "change our genetic")

# dictionary
dict <- list(
  old_age = c("elderly", "age", "over_65", "senior_citizen", "elder", "elders", "old", "over_70"),
  life_expectancy = c("young", "life_expectancy"),
  occupation = c("occupation", "job", "doctor", "doctors", "work", "workers", "profession", "nurse", "physician", 
                 "frontline", "career", "teacher", "teachers", "essential", "frontline", "world_leaders", "system_relevant", "police"),
  preconditions = c("preconditions", "precondition", "health", "conditions", "history", "immune", "disabled", "illness", 
                    "pre-illness", "pre-existing", "ill", "at-risk", "comorbidities", "comorbid", "chronic"),
  has_children = c("has_children", "children_to_raise", "had_kids", "have_children", "with_children", "with_kids", 
                   "has_younger_children", "has_small_children", "family", "dependent", "parent", "dependents", "parents", 
                   "mother", "father", "family_man", "having_children", "age_of_children"),
  registration = c("registration", "registered", "sign_up", "first_come", "enrolled"),
  exposure = c("exposure", "contact", "contacts"),
  productivity = c("merit", "contribution", "education", "service", "value_in_society", "benefit_to_the_general_public", 
                   "productivity", "function_for_society", "knowledge", "who_contribute", "expected_gains"),
  citizenship = c("citizen", "us_born", "american", "americans", "citizenship", "nationality", "country", "german", 
                  "germans", "brazilian", "brazilians", "polish", "pole", "italian", "italians"),
  income = c("income", "purchasing_power"),
  interest = c("want_it", "want_the_vaccine", "everyone_who_wants"),
  gender = c("male", "female", "gender")
)

dict

# gathering dfm with dictionary hits function  ------------
dfm_with_dictionary <- function(text_lines, dict, pattern, country) {
  dfm <- tokens(text_lines, remove_punct = TRUE, 
                remove_symbols = TRUE, 
                remove_numbers = TRUE) %>%
    tokens_compound(., pattern = phrase(phrases)) %>%
    tokens_remove(stopwords("en")) %>%
    dfm(dictionary = dictionary(dict)) 
  
  dfm_df <- convert(dfm_weight(dfm, scheme = "boolean"), 
                    to = "data.frame") %>% 
    dplyr::mutate(country = country)
  return(dfm_df)
}


## df extraction ----

# national dfs
en_dfm_df <- dfm_with_dictionary(text_lines = en, dict = dict, pattern = phrases, country = "USA")
de_dfm_df <- dfm_with_dictionary(text_lines = de, dict = dict, pattern = phrases, country = "DEU")
br_dfm_df <- dfm_with_dictionary(text_lines = br, dict = dict, pattern = phrases, country = "BRA")
it_dfm_df <- dfm_with_dictionary(text_lines = it, dict = dict, pattern = phrases, country = "ITA")
pl_dfm_df <- dfm_with_dictionary(text_lines = pl, dict = dict, pattern = phrases, country = "POL")

# international df
int_dfm_df <- dplyr::bind_rows(en_dfm_df, de_dfm_df, br_dfm_df, it_dfm_df, pl_dfm_df)

## analysis ------

# extracting values for table

perc_country_df <- int_dfm_df %>% dplyr::group_by(country) %>% 
  summarize_if(is.numeric, mean) %>%
  dplyr::mutate_if(is.numeric, function(x) (x*100)) %>%
  dplyr::mutate_if(is.numeric, round, 1) %>%
  dplyr::mutate_if(is.numeric, function(x) {paste0("(",x,"%)")}) %>%
  t() %>% dplyr::as_tibble(., rownames = "type") %>%
  janitor::row_to_names(row_number = 1)

count_country_df <- int_dfm_df %>% dplyr::group_by(country) %>% 
  summarize_if(is.numeric, sum) %>%
  dplyr::mutate_if(is.numeric, round, 2) %>%
  t() %>% dplyr::as_tibble(., rownames = "type") %>%
  janitor::row_to_names(row_number = 1)

perc_cmbd_df <- int_dfm_df %>%
  summarize_if(is.numeric, mean) %>%
  dplyr::mutate_if(is.numeric, function(x) (x*100)) %>%
  dplyr::mutate_if(is.numeric, round, 1) %>%
  dplyr::mutate_if(is.numeric, function(x) {paste0("(",x,"%)")}) %>%
  t() %>% dplyr::as_tibble(., rownames = "country") %>%
  dplyr::rename("Combined" = "V1")

count_cmbd_df <- int_dfm_df %>% 
  summarize_if(is.numeric, sum) %>%
  dplyr::mutate_if(is.numeric, round, 2) %>%
  t() %>% dplyr::as_tibble(., rownames = "country") %>%
  dplyr::rename("Combined" = "V1")

# combined dfs
count_cmbd_df <- dplyr::left_join(count_cmbd_df, count_country_df) %>%
  mutate_at(c("BRA", "DEU", "ITA", "POL", "USA"), function(x)(as.numeric(x)))

perc_cmbd_df <- dplyr::left_join(perc_cmbd_df, perc_country_df) 

# table
totals <- tibble::tibble(country = "totals", 
                         Combined = nrow(int_dfm_df), 
                         BRA = nrow(br_dfm_df),
                         DEU = nrow(de_dfm_df),
                         ITA = nrow(it_dfm_df),
                         POL = nrow(pl_dfm_df), 
                         USA = nrow(en_dfm_df)) %>% dplyr::mutate_if(is.numeric, as.character)

nr <- tibble::tibble(country = "nr",
                     BRA = sum(br == ""),
                     DEU = sum(de == ""),
                     ITA = sum(it == ""),
                     POL = sum(pl == ""), 
                     USA = sum(en == "")) %>%
  dplyr::mutate(Combined = sum(BRA, DEU, ITA, POL, USA)) %>%
  dplyr::mutate_if(is.numeric, as.character) %>%
  dplyr::select(country, Combined, everything())

supp_table_2_df <- count_cmbd_df %>% 
  dplyr::mutate(order = Combined,
                Combined = paste0(Combined," ", perc_cmbd_df$Combined),
                BRA = paste0(BRA," ", perc_cmbd_df$BRA),
                DEU = paste0(DEU," ", perc_cmbd_df$DEU),
                ITA = paste0(ITA," ", perc_cmbd_df$ITA),
                POL = paste0(POL," ", perc_cmbd_df$POL),
                USA = paste0(USA," ", perc_cmbd_df$USA)) %>%
  dplyr::arrange(desc(order)) %>%
  dplyr::select(-order) %>% 
  dplyr::bind_rows(., nr, totals)

supp_table_2_df

## -----------------------------------------------------------------
## SF8. Conjoint (AMCE) --------------------------------------------
## -----------------------------------------------------------------

# estimation function
est_lm <- function(df){
  m <- lm(pat_pref ~ Age + EarlyRegistration + Gender + Children + Job + Precondition + Citizenship, data = df, weights = weight)
  clse_models <-  vcovCL(m, ~ ResponseId)
  d <- data.frame("est" = coef(m),
                  "se" = sqrt(diag(clse_models)),
                  "names" = names(coef(m))) %>%
    filter(names!="(Intercept)")
  return(d)
} 

# by-country estimate
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_lm(.))

# pooled estimate
df_res_pooled <- est_lm(conjoint_df)
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- bind_rows(df_res_pooled, df_res)
supp_fig_8_df <- dplyr::left_join(plot_helper_df, df_res_all) %>%
  dplyr::mutate(est = ifelse(is.na(est), 0, est),
                se = ifelse(is.na(se), 0, se),
                baseline = ifelse(est == 0 & se ==0 & names != "JobTeacher", "Baseline", "AMCE"),
                levl = as.factor(levl),
                country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil")))
supp_fig_8_df$est[supp_fig_8_df$names == "JobTeacher" & supp_fig_8_df$country != "Germany"] <- NA
supp_fig_8_df$se[supp_fig_8_df$names == "JobTeacher" & supp_fig_8_df$country != "Germany"] <- NA

# ticks and labs
xaxis_ticks <- unique(supp_fig_8_df$position)
xaxis_lab <- unique(paste0(supp_fig_8_df$attr, " | ", supp_fig_8_df$levl)) %>% stringr::str_replace(., "Children", "Has children")
vertical_lines <- c(5,8,11,14,17,24)

supp_fig_8 <- ggplot(supp_fig_8_df,aes(x = position, y = est, ymin = est - 1.96 * se, ymax = est + 1.96 * se, 
                   col = baseline)) + 
  geom_hline(yintercept = 0, alpha = 0.3, col = "black") +
  facet_wrap(~ country) + 
  geom_pointrange(size=0.5) +
  scale_x_continuous(breaks = xaxis_ticks,
                     labels = xaxis_lab) +
  scale_color_manual(values = c("black", "grey")) + 
  geom_vline(xintercept=vertical_lines, col="black", linetype="dashed", size=.4) +
  coord_flip() +  xlab("") + ylab("AMCE") + 
  theme_bw() +
  theme(legend.position = "none", text = element_text(size=14),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_8

ggsave("figures/supp-fig-8-amce-all.png", plot=supp_fig_8, height=8, width=10)


## -----------------------------------------------------------------
## SF9. Subgroup (Institutional trust) -----------------------------
## -----------------------------------------------------------------

# estimation function
est_lmr <- function(d){
  
# model estimates
  m <- estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + Children + Job + Precondition +
                             Citizenship, data = d, weights = weight, clusters = personid)
  
  df_plot <- broom::tidy(m) %>% dplyr::filter(term != "(Intercept)")
  
  return(df_plot)
} 

df_res <- conjoint_df %>% 
  group_by(country, trust_cat_institutions_pca) %>%
  do(est_lmr(.))

# pooled estimate
df_res_pooled <- conjoint_df %>% 
  group_by(trust_cat_institutions_pca) %>%
  do(est_lmr(.))
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- bind_rows(df_res_pooled, df_res)

df_names <- tidyr::crossing(plot_helper_df %>% dplyr::distinct(names, attr, levl), df_res_all %>% dplyr::ungroup() %>% tidyr::expand(country, trust_cat_institutions_pca))

figsf9_df <- dplyr::left_join(df_names, df_res_all, by = c("names" = "term", "country", "trust_cat_institutions_pca")) %>% 
  dplyr::mutate(est = ifelse(is.na(estimate), 0, estimate),
                se = ifelse(is.na(std.error), 0, std.error),
                baseline = ifelse(est == 0 & se ==0 & names != "JobTeacher", "Baseline", "ACME"),
                levl = as.factor(levl),
                country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil"))) %>%
  dplyr::filter(!is.na(trust_cat_institutions_pca)) %>%
  dplyr::left_join(plot_helper_df,.)

supp_fig_9 <- ggplot(figsf9_df,aes(x = position)) + 
  geom_hline(yintercept = 0, alpha = 0.3, col = "black") +
  facet_grid(~ factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil"))) + 
  geom_point(data = dplyr::filter(figsf9_df,baseline == "Baseline"), aes(y = 0), shape = 21, color = "black", fill= "white" ,size = 1.5) +
  geom_pointrange(data = dplyr::filter(figsf9_df,baseline == "ACME"), aes(y = est, ymin = est - 1.96 * se, ymax = est + 1.96 * se, color = trust_cat_institutions_pca, shape = trust_cat_institutions_pca), size = 0.3, position = position_dodge(width = 1)) +
  scale_x_continuous(breaks = xaxis_ticks,
                     labels = xaxis_lab) + 
  scale_color_manual(values = c("#bdbdbd", "#737373", "#000000")) +
  geom_vline(xintercept=vertical_lines, col="black", linetype="dashed", size=.4) +
  coord_flip() +  labs(x = "", y = "AMCE")+ 
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank(),
        text = element_text(size=14),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text = element_text(color = "darkblue", face = "bold"
        ),
        axis.text.x = element_text(size = 7))

supp_fig_9

ggsave("figures/supp-fig-9-subgroup-trust.png", plot=supp_fig_9, height=8, width=10)


## -----------------------------------------------------------------
## SF10. Subgroup (Political ideology) -----------------------------
## -----------------------------------------------------------------

df_res <- conjoint_df %>% 
  group_by(country, political_lean) %>%
  do(est_lmr(.))

# pooled estimate
df_res_pooled <- conjoint_df %>% 
  group_by(political_lean) %>%
  do(est_lmr(.))
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- bind_rows(df_res_pooled, df_res)

df_names <- tidyr::crossing(plot_helper_df %>% dplyr::distinct(names, attr, levl), df_res_all %>% dplyr::ungroup() %>% tidyr::expand(country, political_lean))

figsf10_df <- dplyr::left_join(df_names, df_res_all, by = c("names" = "term", "country", "political_lean")) %>% 
  dplyr::mutate(est = ifelse(is.na(estimate), 0, estimate),
                se = ifelse(is.na(std.error), 0, std.error),
                baseline = ifelse(est == 0 & se ==0 & names != "JobTeacher", "Baseline", "ACME"),
                levl = as.factor(levl),
                country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil"))) %>%
  dplyr::filter(!is.na(political_lean)) %>%
  dplyr::left_join(plot_helper_df,.)


supp_fig_10 <- ggplot(figsf10_df,aes(x = position)) + 
  geom_hline(yintercept = 0, alpha = 0.3, col = "black") +
  facet_grid(~ factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil"))) + 
  geom_point(data = dplyr::filter(figsf10_df,baseline == "Baseline"), aes(y = 0), shape = 21, color = "black", fill= "white" ,size = 1.5) +
  geom_pointrange(data = dplyr::filter(figsf10_df,baseline == "ACME"), aes(y = est, ymin = est - 1.96 * se, ymax = est + 1.96 * se, color = political_lean, shape = political_lean), size = 0.3, position = position_dodge(width = 1)) +
  scale_x_continuous(breaks = xaxis_ticks,
                     labels = xaxis_lab) + 
  scale_color_manual(values = c("#bdbdbd", "#737373", "#000000")) +
  geom_vline(xintercept=vertical_lines, col="black", linetype="dashed", size=.4) +
  coord_flip() +  labs(x = "", y = "AMCE")+ 
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank(),
        text = element_text(size=14),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_10

ggsave("figures/supp-fig-10-subgroup-political-lean.png", plot=supp_fig_10, height=8, width=10)


## -----------------------------------------------------------------
## SF11. Citizen-political lean (AMCIE) ----------------------------
## -----------------------------------------------------------------

est_lm_cond <- function(d){
  
# model estimates
  m <- estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + Children + Job + Precondition +
                             Citizenship*political_lean, data = d, weights = weight, clusters = personid)
  
# create estimates
  df_plot <- summary(margins::margins(m, variables = "Citizenship", at = list(political_lean = c("Left","Moderate","Right"))))
  
  return(df_plot)
} 

# estimation countrywise
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_lm_cond(.))

# pooled estimate
df_res_pooled <- conjoint_df %>% 
  do(est_lm_cond(.))
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- dplyr::bind_rows(df_res_pooled, df_res)

# Prepare Df
df_plot <- df_res_all %>%
  mutate(country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil")))

# Figure
supp_fig_11 <- ggplot(df_plot) + 
  geom_hline(yintercept = 0,alpha=0.3,col="red", linetype = "dashed") +
  geom_pointrange(aes(x=political_lean,y=AME,ymin=lower,ymax=upper),
                  position = position_dodge(0.4), size=.8) +
  theme_bw() +
  coord_flip() + facet_wrap(~ country) +  theme(text = element_text(size=12)) +
  xlab("") + ylab("AMCIE") +
  theme(legend.position = "none", axis.text.x = element_text(size=7),
        panel.grid.minor.y =  element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_11

ggsave("figures/supp-fig-11-amcie-citizenship-pol.png", plot=supp_fig_11, width = 6, height = 5)

## -----------------------------------------------------------------
## SF12. Precondition-precondition (AMCIE) -------------------------
## -----------------------------------------------------------------

est_lm_cond <- function(d){
  
# model estimates
  m <- estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + Children + Job + Precondition*resp_precondition +
                             Citizenship, data = d, weights = weight, clusters = personid)
  
# create estimates
  df_plot <- summary(margins::margins(m, variables = "Precondition", at = list(resp_precondition = c("Yes","No"))))
  
  return(df_plot)
} 

# Estimation countrywise
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_lm_cond(.))

# pooled estimate
df_res_pooled <- conjoint_df %>% 
  do(est_lm_cond(.))
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- dplyr::bind_rows(df_res_pooled, df_res)

# Prepare Df
df_plot <- df_res_all %>%
  mutate(country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil")))

# Figure
supp_fig_12 <- ggplot(df_plot) + 
  geom_hline(yintercept = 0,alpha=0.3,col="red", linetype = "dashed") +
  geom_pointrange(aes(x=resp_precondition,y=AME,ymin=lower,ymax=upper),
                  position = position_dodge(0.4), size=.8) +
  theme_bw() +
  coord_flip() + facet_wrap(~ country) +  theme(text = element_text(size=12)) +
  xlab("") + ylab("AMCIE") +
  theme(legend.position = "none", axis.text.x = element_text(size=7),
        panel.grid.minor.y =  element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_12

ggsave("figures/supp-fig-12-amcie-precond-precond.png", plot=supp_fig_12, width = 6, height = 5)

## -----------------------------------------------------------------
## SF13. Has children-parent (AMCIE) -------------------------------
## -----------------------------------------------------------------

est_lm_cond <- function(d){
  
  # Model Estimates
  m <- estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + Children*children_household + Job + Precondition +
                             Citizenship, data = d, weights = weight, clusters = personid)
  
  # Create Estimates
  df_plot <- summary(margins::margins(m, variables = "Children", at = list(children_household = c("Yes","No"))))
  
  return(df_plot)
} 

# Estimation countrywise
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_lm_cond(.))

# pooled estimate
df_res_pooled <- conjoint_df %>% 
  do(est_lm_cond(.))
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- dplyr::bind_rows(df_res_pooled, df_res)

# Prepare Df
df_plot <- df_res_all %>%
  mutate(country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil")))

# Figure
supp_fig_13 <- ggplot(df_plot) + 
  geom_hline(yintercept = 0,alpha=0.3,col="red", linetype = "dashed") +
  geom_pointrange(aes(x=children_household,y=AME,ymin=lower,ymax=upper),
                  position = position_dodge(0.4), size=.8) +
  theme_bw() +
  coord_flip() + facet_wrap(~ country) +  theme(text = element_text(size=12)) +
  xlab("") + ylab("AMCIE") +
  theme(legend.position = "none", axis.text.x = element_text(size=7),
        panel.grid.minor.y =  element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_13

ggsave("figures/supp-fig-13-amcie-child-parent.png", plot=supp_fig_13, width = 6, height = 5)

## -----------------------------------------------------------------
## SF14. Gender-gender (AMCIE) -------------------------------------
## -----------------------------------------------------------------

est_lm_cond <- function(d){
  
  # Model Estimates
  m <- estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender*resp_gender + Children + Job + Precondition +
                             Citizenship, data = d, weights = weight, clusters = personid)
  
  # Create Estimates
  df_plot <- summary(margins::margins(m, variables = "Gender", at = list(resp_gender = c("Male","Female"))))
  
  return(df_plot)
} 

# Estimation countrywise
df_res <- conjoint_df %>% 
  group_by(country) %>%
  do(est_lm_cond(.))

# pooled estimate
df_res_pooled <- conjoint_df %>% 
  do(est_lm_cond(.))
df_res_pooled$country <- "Pooled"

# join results
df_res_all <- dplyr::bind_rows(df_res_pooled, df_res)

# Prepare Df
df_plot <- df_res_all %>%
  mutate(country = factor(country, levels = c("Pooled","USA","Germany","Italy","Poland","Brazil")))

# Figure
supp_fig_14 <- ggplot(df_plot) + 
  geom_hline(yintercept = 0,alpha=0.3,col="red", linetype = "dashed") +
  geom_pointrange(aes(x=resp_gender,y=AME,ymin=lower,ymax=upper),
                  position = position_dodge(0.4), size=.8) +
  theme_bw() +
  coord_flip() + facet_wrap(~ country) +  theme(text = element_text(size=12)) +
  xlab("") + ylab("AMCIE") +
  theme(legend.position = "none", axis.text.x = element_text(size=7),
        panel.grid.minor.y =  element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"
        ))

supp_fig_14

ggsave("figures/supp-fig-14-amcie-gender-gender.png", plot=supp_fig_14, width = 6, height = 5)

## -----------------------------------------------------------------
## ST4. Regression output conjoint ---------------------------------
## -----------------------------------------------------------------

#regression table
models <- list(
  "Pooled" = estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + `Has children` + Job + Precondition + Citizenship, data = conjoint_df, weights = weight, clusters = personid),
  "USA"    = estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + `Has children` + Job + Precondition + Citizenship, data = dplyr::filter(conjoint_df,Q_Language == "EN"), weights = weight, clusters = personid),
  "Germany"= estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + `Has children` + Job + Precondition + Citizenship, data = dplyr::filter(conjoint_df,Q_Language == "DE"), weights = weight, clusters = personid),
  "Italy"  = estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + `Has children` + Job + Precondition + Citizenship, data = dplyr::filter(conjoint_df,Q_Language == "IT"), weights = weight, clusters = personid),
  "Poland" = estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + `Has children` + Job + Precondition + Citizenship, data = dplyr::filter(conjoint_df,Q_Language == "PL"), weights = weight, clusters = personid),
  "Brazil" = estimatr::lm_robust(pat_pref ~ Age + EarlyRegistration + Gender + `Has children` + Job + Precondition + Citizenship, data = dplyr::filter(conjoint_df,Q_Language == "PT-BR"), weights = weight, clusters = personid)
)

modelsummary::modelsummary(models, fmt = 2,
                           statistic = 'conf.int')

texreg::texreg(models, file = "figures/supp-table-4-regression-output.tex", include.rmse = F, fontsize = "footnotesize", ci.test = NA)

## -----------------------------------------------------------------
## SF15. Ranking bars (Pooled) -------------------------------------
## -----------------------------------------------------------------

supp_fig_15 <- ranking_df %>%
  dplyr::mutate(type = factor(type, levels = c("Workers in the\nhealthcare system","People with medical\npreconditions","The elderly","People who have a lot of contact\nwith other people at work",
                                               "Children", "Business leaders", "Leading politicians", "Professional athletes"))) %>%
  dplyr::filter(country == "Pooled" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(x = "Ranking") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%"),
                color = !(type %in% c("Workers in the\nhealthcare system","People with medical\npreconditions","The elderly"))),
            position=position_stack(vjust = 0.5)) +
  theme_bw() +
  scale_y_continuous(labels = scales::label_percent(scale=1)) +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  scale_color_manual(guide = "none", values = c("white", "black")) +
  theme(legend.position = "bottom", legend.title = element_blank(),axis.title.y = element_blank()) + guides(fill = guide_legend(nrow = 1))

supp_fig_15

ggsave("figures/supp-fig-15-ranking-bar-pooled.png", plot=supp_fig_15, width = 12, height = 8)

## -----------------------------------------------------------------
## SF16. Ranking bars (Country-wise) -------------------------------
## -----------------------------------------------------------------

## Pooled 
a <- ranking_df %>%
  dplyr::filter(country == "Pooled" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(y = "Percentage of Participants",
       x = "Ranking",
       title = "(a) Pooled\n") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%")),
            position=position_dodge(width = 0.1), vjust=-0.25, 
            color = "black") +
  theme_bw() +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  facet_wrap(~type) +
  theme(legend.position = "none", text = element_text(size=12),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"),
        plot.title = element_text(color = "darkblue", face = "bold", hjust = 0.5)
        
  ) +
  scale_y_continuous(limits = c(0,70))

## United States
b <- ranking_df %>%
  dplyr::filter(country == "USA" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(y = "Percentage of Participants",
       x = "Ranking",
       title = "(b) USA\n") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%")),
            position=position_dodge(width = 0.1), vjust=-0.25, 
            color = "black") +
  theme_bw() +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  facet_wrap(~type) +
  theme(legend.position = "none", text = element_text(size=12),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"),
        plot.title = element_text(color = "darkblue", face = "bold", hjust = 0.5)
        
  ) +
  scale_y_continuous(limits = c(0,70))

## Germany
c <- ranking_df %>%
  dplyr::filter(country == "Germany" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(y = "Percentage of Participants",
       x = "Ranking",
       title = "(c) Germany\n") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%")),
            position=position_dodge(width = 0.1), vjust=-0.25, 
            color = "black") +
  theme_bw() +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  facet_wrap(~type) +
  theme(legend.position = "none", text = element_text(size=12),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"),
        plot.title = element_text(color = "darkblue", face = "bold", hjust = 0.5)
        
  ) +
  scale_y_continuous(limits = c(0,70))

## Italy
d <- ranking_df %>%
  dplyr::filter(country == "Italy" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(y = "Percentage of Participants",
       x = "Ranking",
       title = "(d) Italy\n") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%")),
            position=position_dodge(width = 0.1), vjust=-0.25, 
            color = "black") +
  theme_bw() +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  facet_wrap(~type) +
  theme(legend.position = "none", text = element_text(size=12),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"),
        plot.title = element_text(color = "darkblue", face = "bold", hjust = 0.5)
        
  ) +
  scale_y_continuous(limits = c(0,70))

## Poland
e <- ranking_df %>%
  dplyr::filter(country == "Poland" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(y = "Percentage of Participants",
       x = "Ranking",
       title = "(e) Poland\n") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%")),
            position=position_dodge(width = 0.1), vjust=-0.25, 
            color = "black") +
  theme_bw() +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  facet_wrap(~type) +
  theme(legend.position = "none", text = element_text(size=12),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"),
        plot.title = element_text(color = "darkblue", face = "bold", hjust = 0.5)
        
  ) +
  scale_y_continuous(limits = c(0,70))

## Brazil
f <- ranking_df %>%
  dplyr::filter(country == "Brazil" & !is.na(val)) %>%
  dplyr::count(type, val) %>%       
  dplyr::group_by(type) %>%
  dplyr::mutate(pct= prop.table(n) * 100) %>%
  ggplot() + aes(factor(val), pct, fill=type) +
  geom_bar(stat="identity", color = "black") +
  labs(y = "Percentage of Participants",
       x = "Ranking",
       title = "(f) Brazil\n") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct),"%")),
            position=position_dodge(width = 0.1), vjust=-0.25, 
            color = "black") +
  theme_bw() +
  scale_fill_manual(values = rev(c('#b35806','#e08214','#fdb863','#fee0b6','#d8daeb','#b2abd2','#8073ac','#542788'))) +
  facet_wrap(~type) +
  theme(legend.position = "none", text = element_text(size=12),
        panel.grid.minor.y =  element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(color = "darkblue", face = "bold"),
        plot.title = element_text(color = "darkblue", face = "bold", hjust = 0.5)
        
  ) +
  scale_y_continuous(limits = c(0,70))

#pooled
ggsave("figures/supp-fig-16-a-ranking-bar-pooled.png", plot=a, width = 12, height = 8)

#usa
ggsave("figures/supp-fig-16-b-ranking-bar-usa.png", plot=b, width = 12, height = 8)

#ger
ggsave("figures/supp-fig-16-c-ranking-bar-ger.png", plot=c, width = 12, height = 8)

#italy
ggsave("figures/supp-fig-16-d-ranking-bar-ita.png", plot=d, width = 12, height = 8)

#poland
ggsave("figures/supp-fig-16-e-ranking-bar-pol.png", plot=e, width = 12, height = 8)

#brazil
ggsave("figures/supp-fig-16-f-ranking-bar-bra.png", plot=f, width = 12, height = 8)
